home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / bandred.r < prev    next >
Encoding:
Text File  |  1994-02-21  |  1.7 KB  |  70 lines  |  [TEXT/RLAB]
  1. //
  2. // BANDRED  B = BANDRED(A, KL, KU) is a matrix orthogonally equivalent to A
  3. //          with lower bandwidth KL and upper bandwidth KU
  4. //          (i.e. B(i,j) = 0 if i>j+KL or j>i+KU).
  5. //          The reduction is performed using Householder transformations.
  6. //          If KU is omitted it defaults to KL.
  7.  
  8. //          This routine is called by RANDSVD.
  9. //          This is a `standard' reduction.  Cf. reduction to bidiagonal form
  10. //          prior to computing the SVD.  This code is a little wasteful in that
  11. //          it computes certain elements which are immediately set to zero!
  12.  
  13. rfile house
  14.  
  15. bandred = function(A, kl, ku)
  16. {
  17.   local(a, beta, flip, j, ltmp, m, n, temp, v, z);
  18.  
  19.   if(!exist (ku)) { ku = kl; }
  20.  
  21.   if( kl == 0 && ku == 0 ) {
  22.     error("You've asked for a diagonal matrix. In that case use the SVD!");
  23.   }
  24.  
  25.   // Check for special case where order of left/right transformations matters.
  26.   // Easiest approach is to work on the transpose, flipping back at the end.
  27.  
  28.   a = A;
  29.   flip = 0;
  30.   if( ku == 0 )
  31.   {
  32.     a = a';
  33.     temp = kl; 
  34.     kl = ku;
  35.     ku = temp; 
  36.     flip = 1;
  37.   }
  38.  
  39.   m = size(a)[1]; n = size(a)[2];
  40.  
  41.   if( (z=min( [min([m,n]),max([m-kl-1,n-ku-1]) ] )) >= 1 )
  42.   {
  43.     for( j in 1:z )
  44.     {
  45.       if( j+kl+1 <= m )
  46.       {
  47.          ltmp = house( a[j+kl:m; j] );
  48.          beta = ltmp.beta; v = ltmp.v;
  49.          temp = a[j+kl:m; j:n];
  50.          a[j+kl:m; j:n] = temp - beta*v*(v'*temp);
  51.          a[j+kl+1:m; j] = zeros(m-j-kl,1);
  52.       }
  53.  
  54.       if( j+ku+1 <= n )
  55.       {
  56.          ltmp = house( a[j; j+ku:n]' );
  57.          beta = ltmp.beta; v = ltmp.v;
  58.          temp = a[j:m; j+ku:n];
  59.          a[j:m; j+ku:n] = temp - beta*(temp*v)*v';
  60.          a[j; j+ku+1:n] = zeros(1,n-j-ku);
  61.       }
  62.     }
  63.   }
  64.  
  65.   if( flip ) {
  66.    a = a';
  67.   }
  68.   return a;
  69. };
  70.